#
#  figure7a.r
#  This code creates the first row of four graphs in figure 7, showing the distribution of imputations over
#    missing values, compared to the observed data, for the key variable Welfare spending as a percent of GDP.
#
#  23/2/10, jH
#

stub<-"p3ny"     # this is the stub name of the imputed datasets

nimp<-100        # number of imputed datasets to cycle over
alpha<-0.95      # this drives the size of confidence interval (but is strictly the "alpha" value of the ci)
jitter<-0.2      # a graphical parameter


a<-load("burgoonsubset.rdata")
impdata<-reasonabledata

subcty<-c(53,155,452,395)                               # These are the countries to display
countrynames<-c("Barbados","Chile","Ghana","Iceland")

viewedlimits<-matrix(c(8,18,8,25,0,12,9,23),2,4)        # A series of axis limits
viewedlimits<-t(viewedlimits)


ci.upper<-function(data){
  alpha<-0.95   #local alpha (internal to function), not global alpha, makes apply easier, and yet cumbersome
  a<-sort(data)
  b<-a[round(length(a)*alpha)] 
  return(b)
}

ci.lower<-function(data){
  alpha<-0.95
  a<-sort(data)
  b<-a[round(length(a)*(1-alpha))] 
  return(b)
}

period<-c(1974,1996)

yearseq<-(period[1]+1):(period[2]-1)
xyears<-rbind(yearseq,yearseq,NA)
delta<-0.1
xyearsdelta<-rbind(yearseq-delta,yearseq+delta,NA)

l<-1
w<-4
lw<-l*w
par(mfrow=c(l,w)) 




for(i in 1:length(subcty)){          # Loop over countries

  flag<-impdata$cow==subcty[i] & impdata$year>period[1] & impdata$year<period[2]
  tempdata<-impdata[flag, ]
  allmain<-countrynames[i]
  allylim<-viewedlimits[i,]
 
  r1<-is.na(tempdata$welfarepercent)
  stackdata<-matrix(0,period[2]-period[1]-1,nimp)
  impeddata<-matrix(0,period[2]-period[1]-1,nimp)

  for(j in 1:nimp){
    temp.imp<-read.csv(file=paste(stub,j,".csv",sep=""))
    impflag<-temp.imp$cow==subcty[i] & temp.imp$year>period[1] & temp.imp$year<period[2]
    impeddata[,j]<-temp.imp$welfarepercent[impflag]
  }

  for(i in 1:nrow(stackdata)){
    stackdata[i,]<-sort(impeddata[i,])
  }

#  s.upper <-apply(stackdata,1,ci.upper)
#  s.median<-apply(stackdata,1,median)
#  s.lower <-apply(stackdata,1,ci.lower)

   cilr<-c(round(nimp*(1-alpha)),round(nimp*alpha))

   missing<-is.na(tempdata$welfarepercent)
   years.missing<-xyears[1,missing]


   s.upper <-  stackdata[,cilr[1]]
   s.median <- stackdata[,round(nimp/2)]
   s.lower <-  stackdata[,cilr[2]]

  # Confidence Intervals  
  plot(xyears[,r1],rbind(s.upper[r1],s.lower[r1],NA),main=allmain,xlab="Year",ylab="Welfare Percent",ylim=allylim,xlim=period,type="l",col="slategrey")
  par(new=TRUE)
  plot(xyearsdelta[,r1],rbind(s.median[r1],s.median[r1],NA),main=allmain,xlab="Year",ylab="Welfare Percent",ylim=allylim,xlim=period,type="l",col="slategrey")
  par(new=TRUE)

  # Some (m) Imputed Values
  m<-5
  plot((years.missing*matrix(1,length(years.missing),m)),impeddata[missing,1:m],main=allmain,xlab="Year",ylab="Welfare Percent",ylim=allylim,xlim=period,type="p",col="blue")
  par(new=TRUE)

  # Observed Data
  plot(tempdata$year,tempdata$welfarepercent,main=allmain,xlab="Year",ylab="Welfare Percent",ylim=allylim,xlim=period,pch=16)
  
  if(round(i/lw)*lw==i){
    par(ask=TRUE)
  }else{
    par(ask=FALSE)
  }

}

  
